---
title: "To What Extent do Political Parties Influence Public Opinion?"
author: "Evelyn Cai, Genevieve Waldorf, and Yao Yu"
date: "12/14/2020"
output:
  pdf_document: default
always_allow_html: true
bibliography: bib.bib
nocite: '@*'
fontsize: 12pt
indent: true
header-includes:
    - \usepackage{setspace}\doublespacing
abstract: "In this paper, we examine the extent to which political parties influence public opinion for voters that identify with the party. We examine Danish panel survey data from 2010-2011 collected by Slothuus and Bisgaard (2020) to determine the impact of party policy changes towards unemployment benefits and early retirement on public support. We employ the coarsened exact matching technique as a preprocessing tool to reduce bias on matched control and treatment units. We find that party-identifying citizens’ support levels for certain policies are influenced by the official political party’s support to a lesser extent than previously shown, especially for the retirement policy. We also consider evidence for reverse causation."
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)

##########################################
# HOW POLITICAL PARTIES SHAPE PUBLIC OPINION IN THE REAL WORLD
# AUTHORS: Martin Bisgaard & Rune Slothuus
# VERSION: SEPTEMBER 2020 
# THIS SCRIPT: REPRODUCES ANALYSES IN MAIN ARTICLE
###########################################

rm(list=ls())

#Set user specific directory
#setwd("/mydirectory")

#Sys.setenv(LANGUAGE="en")
#Sys.setlocale("LC_TIME", "English") 


#installing required packages
#install.packages("rio")
#install.packages("foreign")
#install.packages("car")
#install.packages("plyr")
#install.packages("plm")
#install.packages("lme4")
#install.packages("lmtest")
#install.packages("stargazer")
#install.packages("xtable")
#install.packages("psych")
#install.packages("grDevices")
#install.packages("plyr")
#install.packages("reshape2")
#install.packages("pBrackets")
#install.packages("MatchIt")
#install.packages("sandwich")
#install.packages("knitr)
#install.packages("kableExtra)


#loading required packages
require(rio)
require(foreign)
require(car)
require(plyr)
require(plm)
require(lme4)
require(lmtest)
require(stargazer)
require(xtable)
require(psych)
require(grDevices)
require(plyr)
require(reshape2)
require(pBrackets)
require(MatchIt)
require(sandwich)


#load panel data
dat.balance<-read.dta("paneldata_long.dta",convert.factors=TRUE)
paneldat<-read.dta("paneldata_wide.dta",convert.factors=TRUE)
dat.balance$partyAll<-recode(dat.balance$partyAll,'"NA"=NA') 
#fixing NA value on party variable 

##loaded dataset is unbalanced. Store it as separate object to load it easily
##for future use.
dat<-dat.balance
ogdat <- dat
```

# Introduction

Party politics play out in every type of regime: from competitive authoritarian regimes that mimic democratic processes to actual democratic states. It is of substantive interest to examine the influence of political parties on public opinion, as it provides a more nuanced understanding of the relationship between political elites and their constituents. We know that constituents generally vote ideologically, or how well they think a politician will do to pass policy that the constituent views favorably (Degan & Merlo, 2009). Can politicians influence what constituents deem favorable? If so, the implications are large for a better understanding of individual politicians, party organizations, and democratic processes in general.

In *How Political Parties Shape Public Opinion in the Real World*, Slothuus and Bisgaard employ a panel study to examine the impact of Danish parties on the policy support of their constituents. In negotiations with other governing parties such as the Danish People’s Party (DPP), the liberal Venstre party proposed to cut unemployment benefits and the early retirement program in response to the deepening economic crisis in May 2010 and January 2011 respectively. Through analyzing the levels of policy support before and after these party-proposed changes in treatment groups (those who identified with the DPP or Venstre) and control groups (those who identified with other parties that did not undergo a platform change) using propensity score matching, Slothuus and Bisgaard find that policy support substantially increased by about 15 and 17 percentage points for each respective policy. They found about the same increase for each policy, which provides evidence that these Danish parties did succeed in influencing their constituents.

However, propensity score matching (PSM) generates imbalance and thus results in data points that make causal inference difficult. As more observations are pruned to make better “matches” with more propensity-score balance, PSM results in more model dependence and general imbalance because PSM attempts to squash k-dimensions into one using a logit link, which allows for previously imbalanced covariates to contribute to further imbalance. This imbalance is a consequence that potentially contributes to more bias and model dependence.

Therefore, we will examine the quality of matches generated through coarsened exact matching (CEM), which attempts to nonparametrically match treatment and control units 1-to-1. CEM uses adjustable bins to temporarily coarsen each independent variable and prunes away observations with these coarsened bins, therefore enabling increased capture of information with more units than exact matching. However, there are distance trade-offs in CEM.

Similar to setting a caliper in PSM, the strictness of cut-points are customizable. There is a bias-variance tradeoff in matching, which was especially apparent in our application of coarsened exact matching. As the number of pruned units increases, bias decreases, but variance increases. As such, the confidence intervals are wider for CEM-matched units than PSM-matched units. Despite the overlap of some confidence intervals for the CEM-matched units and the PSM-matched units with this overlap between propensity and coarsened matching clearly describes substantive difference between the two methods.


# Our Extension

Bisgaard and Slothuus (2020) present five figures and two tables in their main paper. Our contribution focuses primarily on their figures 3 and 4, which show the change in public opinion among survey respondents after the Liberal Party and Danish People’s Party (DPP) Party shift their positions on two policies. Figure 3 in Bisgaard and Slothuus (2020) examines the changes in public opinion on the policy shift to decrease unemployment benefits from four to two and a half years. Figure 4 in Bisgaard and Slothuus (2020) examines the changes in the policy shift to abolish early retirement. 

In our CEM matching process, we kept the same variables that Bisgaard and Slothuus (2020) used to match in their PSM process. Our defined cut-points remained constant for both sets of data relating to the policy changes. We chose specific cut-points that would place respondents into natural groups, defined below. The cut points made for each covariate were based on substantive research of each variable and analysis of their distribution among treated and non-treated groups. 

Pressure is an index of variables from 0 to 1 used to measure whether the respondent thought that the welfare state was under economic pressure. In the data, we found that the actual values ranged from 0 to 0.8, so we placed our cut-point at the half-way point of 0.4. The variable that asked for the respondents’ age ranged from 0 to 100, so we placed cut-points at 20, 40, 60, and 80. The variable that asked for the respondents’ income ranged from 1 to 9 and increments by DKK 100,000 with 9 as “Don't know / refuse to answer”. We found the median income response was 6, so we used that as our cut-point. Education responses range from 1 to 7 at different levels of schooling. We chose to group 1 through 4 and 5 through 7, the cut-point between high school and college/ higher education. 


# Unemployment Benefit Period Cut

Tables 1 and 2 provide motivation for the improvement of matching. The first pair of observations, although matched considerably well on other covariates, differ greatly in terms of age, income, benefitW1, and benefitW2. The third pair of matches differs greatly in terms of education, age, income, and benefitW1. Through these observations, it is evident that there is also inconsistency regarding which covariate is matched inaccurately. These matching discrepancies demonstrate the shortcomings of this method of matching. In contrast, coarsened exact matching creates pairs of observations whose difference between covariates are consistently smaller. The pairs of observations have exact matches on covariates of education, sex, unemployed, income, benefitW1, and benefitW2. For the covariates in which there was not an exact match, the difference between the pairs of observations were consistently smaller than the differences present in the PSM observations. Through substantive research and analysis of covariate distribution, we were able to establish cut points in the coarsened exact matching process that produced these results. 



```{r, message=FALSE}
# message=FALSE to supress tidyverse conflicts

library(MatchIt)


########################
# Reproducing Figure 3 #
########################

#balancing panel
dat.balance<-dat
dat.balance<-dat.balance[!dat.balance$id %in% (dat.balance[is.na(dat.balance$benefitClean)|is.na(dat.balance$partyAll),]$id),]
# dim(dat.balance)

## matched control group ##
	covariates<-c("pressure","education","sex","unemployed","age","income",
	              "benefitW1","benefitW2")
	dat_match<-dat.balance[,c(covariates,"partyAll","id","time")]
	dat_match<-dat_match[unique(dat_match$id)&dat_match$time==1,] #return to "wide" format
	dat_match<-na.omit(dat_match)
	dat_match$treated<-recode(dat_match$partyAll,'"V"=1;"DF"=1;else=0') 
	
	#code treatment group vs. control pool
#	table(dat_match$treated)

	#get rid of danish letters on values on covariates
	dat_match$income<-as.numeric(as.factor(dat_match$income))
	#dat_match$education<-as.factor(as.numeric(as.factor(dat_match$education)))
	dat_match$education<-(as.numeric(as.factor(dat_match$education)))

	match <- matchit(treated ~factor(education)+pressure+unemployed+sex+age+factor(income)+
	                   factor(benefitW1)+factor(benefitW2),
                     method = "nearest", data = dat_match)

	press <- c(0, 0.4, 0.8)
	
	ag <- c(0,20,40,60,80,100)
	
	inc <- c(1,6,9)

	educ.grp <- list(c("1", "2", "3", "4"), c("5", "6", "7"))
	
	ret.grp <- list(c("1", "2"), c("3","4","5"))

	inc.grp <- list(c("1", "2", "4", "5"), c("6","7","8", "9"))
	
	s <- list(c("0"), c("1"))
	
	emp <- list(c("0"), c("1"))
	
	educ <- c(1,5,7)
	
	ret <- c(1,3,5)
	
	match2 <- matchit(treated ~factor(education)+pressure+unemployed+sex+age+factor(income)+
	                   factor(benefitW1)+factor(benefitW2),data = dat_match,
                    method = "cem", cutpoints = list(age = ag, 
                    pressure = press), 
                    grouping = list(c(benefitW1 = ret.grp, benefitW2 = ret.grp,
                    education = educ.grp, income = inc.grp, sex = s, unemployed = emp)), k2k = TRUE)
	
	
	
	

	dta_m<-match.data(match)
	matchID<-dta_m$id[dta_m$treated==0] #store ID of matched control group

	dta_m2<-match.data(match2)
	matchID2<-dta_m2$id[dta_m2$treated==0] #store ID of matched control group

## calculate quantities of interest for plotting

dat.balance$partygroup<-recode(dat.balance$partyAll,'"V"=1;"DF"=1;else=0')
dat.balance$partygroup<-factor(dat.balance$partygroup)
# table(dat.balance$partygroup)/5


dat.balance$partygroup_matched<-recode(dat.balance$partyAll,'"V"=1;"DF"=1;else=NA')
dat.balance$partygroup_matched[dat.balance$id%in%matchID] <- 0

dat.balance$partygroup_matched2<-recode(dat.balance$partyAll,'"V"=1;"DF"=1;else=NA')
dat.balance$partygroup_matched2[dat.balance$id%in%matchID2] <- 0

#calculate quantities of interest
SE=function(x)sd(x)/sqrt(length(x)-1) #function to calculate naive SEs
out<-aggregate(benefitClean~partygroup*time,data=dat.balance,mean)
out.se<-aggregate(benefitClean~partygroup*time,data=dat.balance,SE)

outm<-aggregate(benefitClean~partygroup_matched*time,data=dat.balance,mean)
outm.se<-aggregate(benefitClean~partygroup_matched*time,data=dat.balance,SE)

outm2<-aggregate(benefitClean~partygroup_matched2*time,data=dat.balance,mean)
outm2.se<-aggregate(benefitClean~partygroup_matched2*time,data=dat.balance,SE)

######
# Showing matches in tables

library(tidyverse)

benefit_psm <- match.data(object = match) %>%
  filter(subclass %in% c("24", "134", "227")) %>%
  select(-time, -treated, -weights) %>%
  mutate(pressure = round(pressure, digits = 2),
         distance = round(distance, digits = 2)) %>%
  rename(group = subclass,
         party = partyAll,
         educ = education,
         unemp = unemployed,
         dist = distance,
         pres = pressure,
         inc = income,
         benW1 = benefitW1,
         benW2 = benefitW2) %>%
  arrange(group) %>% 
  select(group, id, dist, pres, educ, sex, unemp, age, inc, benW1, benW2, party)

table_1 <- knitr::kable(benefit_psm, caption = "PSM Match on Unemployment Policy Responses")

benefit_cem <- match.data(object = match2) %>%
  filter(subclass %in% c("1", "17", "25")) %>%
  select(-time, -treated, -weights) %>%
  mutate(pressure = round(pressure, digits = 2)) %>%
  rename(group = subclass,
         party = partyAll,
         educ = education,
         unemp = unemployed,
         pres = pressure,
         inc = income,
         benW1 = benefitW1,
         benW2 = benefitW2) %>%
  arrange(group) %>% 
  select(group, id, pres, educ, sex, unemp, age, inc, benW1, benW2, party)

table_2 <- knitr::kable(benefit_cem, caption = "CEM Match on Unemployment Policy Responses")
	
```

\newpage

```{r}
table_1
```


```{r}
table_2
```

```{r, include=FALSE}

## PROPENSITY MATCH
plot(summary(match))
```


```{r, include=FALSE}

## CEM MATCH

plot(summary(match2))
```


```{r, include=FALSE}

## UNMATCHED IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

cov_means <- dat_match %>%
  group_by(treated) %>%
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var <- dat_match %>%
  group_by(treated) %>%
   dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre0 <- cov_means %>% 
  left_join(cov_var, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Before matching')


ggplot(bal_pre0, aes(x = var, y = SMD)) + 
  geom_bar(stat = 'identity',
           col = 'black') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()


```


```{r, include=FALSE}

## PROPENSITY SCORE IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

#dat_bal <- dat.balance
#dat_bal$treated<-recode(dat_bal$partyAll,'"V"=1;"DF"=1;else=0')


cov_means <- dta_m %>%
  group_by(treated) %>%
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", "income", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var <- dta_m %>%
  group_by(treated) %>%
   dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", "income", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre <- cov_means %>% 
  left_join(cov_var, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Propensity Score')


ggplot(bal_pre, aes(x = var, y = SMD)) + 
  geom_bar(stat = 'identity',
           col = 'black') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()

```


```{r, include=FALSE}

## PROPENSITY VS CEM IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

#dat_bal <- dat.balance
#dat_bal$treated<-recode(dat_bal$partyAll,'"V"=1;"DF"=1;else=0')


cov_means2 <- dta_m2 %>%
  group_by(treated) %>%
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", "income", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var2 <- dta_m2 %>%
  group_by(treated) %>%
   dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", "income", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre2 <- cov_means2 %>% 
  left_join(cov_var2, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'CEM')


total_bal <- left_join(bal_pre,bal_pre2,by = c('var', 'type'))

total_bal <- rbind(bal_pre,bal_pre2)



ggplot(total_bal, aes(x = var, y = SMD, fill = type)) + 
  geom_bar(stat = 'identity',
           col = 'black', position = 'dodge') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()




```

Figure 1 shows a direct comparison of imbalance between the data before matching, after propensity score matching, and after coarsened exact matching by calculating the normalized difference amongst all covariates. This provides a visual representation of the improvements made on imbalance. The normalized difference values that result from CEM are smaller than the values that result from PSM for all covariates other than unemployment and age. Although there is not improvement of balance on every covariate, the improvements made suggest that when using PSM, the analysis could produce biased results. This implies that the trend of policy support could not be as strong as the original paper demonstrated. 

```{r, warning=FALSE, fig.cap = "Comparison of variable imbalance among the respondents' to unemployment policy. CEM matching results in a better match than the PSM matching used in the original paper."}

# Figure 1

## COMPLETE COMPARISON

library(StatMatch)
## Step 1: MAHALANOBIS DISTANCE MATCHING 
## Covariate Data Matrix 

big_X <- dat_match %>% 
  dplyr::select(c( "pressure", "unemployed", "sex", "age",factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Only treated units 

X_treated <- dat_match %>%
  filter(treated == 1) %>%
  dplyr::select(c( "pressure", "unemployed", "sex", "age",factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Only control units 

X_control <- dat_match %>%
  filter(treated == 0) %>%
  dplyr::select(c( "pressure", "unemployed", "sex", "age",factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Calculate Mahalanobis Distance between all observation-pairs 
## using mahalanobis.dist function from the StatMatch package 

M <- mahalanobis.dist(data.x = X_treated, 
                      data.y = X_control,
                      vc = var(big_X))

## For each treated unit, find closest control unit 
## in terms of Mahalanobis Distance 

match_ids <- unlist(apply(M, 1, function(i){which(i == min(i))[1]}))

## Create matched dataset 

matched_data <- rbind(subset(dat_match, treated == 1),
                      subset(dat_match, treated == 0)[match_ids,])

matched_data_cem <- dta_m2

matched_data_prop <- dta_m


## Check balance after matching ! 

cov_means_matched <- matched_data %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched <- matched_data %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post <- cov_means_matched %>% 
  left_join(cov_var_matched, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Mahalanobis')
#####
cov_means_matched2 <- matched_data_cem %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched2 <- matched_data_cem %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post2 <- cov_means_matched2 %>% 
  left_join(cov_var_matched2, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'CEM')

####
cov_means_matched3 <- matched_data_prop %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c( "pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched3 <- matched_data_prop %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("pressure", "unemployed", "sex", "age", factor("benefitW1"), factor("benefitW2"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post3 <- cov_means_matched3 %>% 
  left_join(cov_var_matched3, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'PSM')

## Compare balance before and after matching visually 


ggplot(rbind(bal_pre0, bal_post2, bal_post3), aes(x = var, y = SMD, 
                                     col = type, group = type, shape = type)) +
  geom_point() + 
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw() + 
  geom_hline(yintercept = 0,
             linetype = 'dotted') + 
  theme(legend.position = 'bottom',
        legend.title = element_blank())

```

When replicating the figure from the original paper (Figure 2), we see that the improvement of matching leads to results that corroborate the author’s initial arguments. The trend of public support from the matches of PSM follows a similar trend to the policy support for the matches of CEM. The confidence intervals from the coarsened exact matching are a result of the aforementioned tradeoffs of this method of matching. However, our study of PSM imbalance and further analysis of the data led to the discovery of the inherently biased dataset. The number of unemployed observations is 745 while the number of employed observations is 13765. As seen by Figure 1 there is no value for CEM for the unemployed covariate. Because there are considerably few unemployed observations to match on, CEM resulted in matches of only employed observations. 

```{r, fig.cap = "Respondents' support for the policy to reduce the unemployment benefit period. Vertical dotted line represents the date in which the treatment party shifted their policy support. Confidence Inervals for each point are at the 95% level."}

# Figure 2

######
# PLOTTING

#sets correct spacing between waves
x<-c("2010-02-22","2010-04-05","2010-06-15","2011-01-11","2011-06-25") 
x<-as.Date(x)


c1="grey70"
c3 = "dodgerblue"

#pdf('benefits_plot.pdf',width=6,height=4)
  par(omd=c(0.05,1,0,1),mai=c(1,1,.1,.1),mfrow=c(1,1))
  y=out[out$partygroup==1,3]
  se=out.se[out.se$partygroup==1,3]
  plot(x,y,ylim=c(.2,.8),xlim=c(x[1],as.Date("2011-08-01")),main="",
       xaxt="n",yaxt="n",ylab="",xlab="")
  segments(x,y+1.96*se,x,y-1.96*se)
  lines(x,y,lty=1)
  points(x,y,pch=21,bg="black")
  y=out[out$partygroup==0,3]
  se=out.se[out.se$partygroup==0,3]
  lines(x,y,lty=1,col=c1)
  points(x,y,pch=24,bg=c1,col=c1)
  segments(x,y+1.96*se,x,y-1.96*se,col=c1)
  y=outm[outm$partygroup==0,3]
  se=outm.se[outm.se$partygroup==0,3]
  lines(x,y,lty=1,col=c1)
  points(x,y,pch=21,bg=c1,col=c1)
  segments(x,y+1.96*se,x,y-1.96*se,col=c1)
  y=outm2[outm2$partygroup==0,3]
  se=outm2.se[outm2.se$partygroup==0,3]
  lines(x,y,lty=1,col=c3)
  points(x,y,pch=21,bg=c3,col=c3)
  segments(x,y+1.96*se,x,y-1.96*se,col=c3)

#redraw CI's of treatment group (for overlap)
  y=out[out$partygroup==1,3]
  se=out.se[out.se$partygroup==1,3]
  segments(x,y+1.96*se,x,y-1.96*se)
  points(x,y,pch=21,bg="black")

  axis(2,las=2)
  axis(2,at=seq(0,1,.025),tck=-.015,labels=NA)
  mtext(2,text="Policy Support",line=3)
  axis.Date(1,at=seq(as.Date("2010-03-01"),as.Date("2011-08-01"), by="month"),format="%b",las=1,cex.axis=1)
  axis.Date(1,at=c(as.Date("2010-03-01"),as.Date("2011-01-01")),format="%Y",
            las=1,cex.axis=1,line=1,tick="FALSE")
  abline(v=as.Date("2010-05-21"),lty=2)

  legend('bottomright',pch=c(21,24,21),pt.bg=c("black",c1,c1, c3),
         col=c("black",c1,c1,c3),legend=c("Treated","Control (all)","Matched (PSM)", "Matched (CEM)"),cex=.8,pt.cex=1,bty='n')

#dev.off()
```

The dataset of this study regarding alterations of unemployment benefits involved participants where less than six percent of the participants were unemployed. This policy directly affects the unemployed population by posing to cut unemployment benefits. The original paper concluded from this analysis that the public opinion followed the policy shift of the parties. This implies that after political parties pose to decrease unemployment benefits, the public follows suit. Without equally considering the opinions of unemployed and employed observations, the proposed outcome could be biased in favor of employed participants. It follows intuition that unemployed participants would be more hesitant to support policies that directly and negatively impact their well being. This presence of imbalance suggests that the current analysis of public support of the policy shift of unemployment benefits is skewed towards the employed public beliefs. To lessen bias, there needs to be balance among employed and unemployed respondents. Further analysis can be done regarding this topic. Our suggestions include: bootstrapping to simulate an increased number of unemployed participants or dropping additional employed participants.

# Abolish Early Retirement

In Bisgaard and Slothuus (2020), they find that Liberals were more likely to support abolishing the early retirement policy after the party shifted their position on the issue from not supporting to supporting. However, we can see that the comparison between the CEM and PSM matched pairs in tables 3 and 4 show that the CEM had a much better match than the PSM. The sample of CEM matches show exact matches for most variables with only age differing slightly. One the other hand, PSM in table 3 may have been able to match every treated respondent, the Liberals, to a control. However, we can see that there is imbalance among these matched pairs. Group 22, for example, is a match between a respondent with polar opposite income levels. Here, we can see that the two respondents have completely different responses in retireW1, retireW2, and retireW2, the three surveys they took at different times in 2010 before the Liberal party shifted their policy position. Using CEM, we were able to reduce the imbalance at the expense of having fewer matched pairs.

```{r, echo=FALSE}
library(tidyverse)
library(MatchIt)
########################
# Reproducing Figure 4 #
########################

#balancing panel for the issue
dat.balance<-dat
dat.balance<-dat.balance[!dat.balance$id %in% (dat.balance[is.na(dat.balance$retire_cutClean)|is.na(dat.balance$partyAll),]$id),]
# dim(dat.balance)

## matched control group ##
	covariates<-c("unemployed","pressure","education","sex","age",
	              "income","retireW1","retireW2","retireW3")
	dat_match<-dat.balance[,c(covariates,"partyAll","id","time")]
	dat_match<-dat_match[unique(dat_match$id)&dat_match$time==1,] 
	#return to "wide" format
	dat_match<-na.omit(dat_match)
	
	dat_match <- dat_match %>%
	  mutate(treated = case_when(partyAll == "V" ~ 1,
	                             TRUE ~ 0))
	#dat_match$treated<-recode(dat_match$partyAll,'"V"=1;else=0') 
	#code treatment group vs. control pool
#	table(dat_match$treated)

	#get rid of danish letters on values on covariates
	dat_match$income<-as.numeric(as.factor(dat_match$income))
	#dat_match$education<-as.factor(as.numeric(as.factor(dat_match$education)))
	dat_match$education<-(as.numeric(as.factor(dat_match$education)))

	#388 matched
	match <- matchit(treated ~ factor(retireW1)+factor(retireW2)+factor(retireW3)+
	                   pressure+sex+unemployed+factor(income)+factor(education)+age,
                     method = "nearest", data = dat_match)
	
	press <- c(0, 0.4, 0.8)
	
	ag <- c(0,20,40,60,80,100)
	
	inc <- c(1,6,9)

	educ.grp <- list(c("1", "2", "3", "4"), c("5", "6", "7"))
	
	ret.grp <- list(c("1", "2"), c("3","4","5"))

	inc.grp <- list(c("1", "2", "4", "5"), c("6","7","8", "9"))
	
	s <- list(c("0"), c("1"))
	
	emp <- list(c("0"), c("1"))
	
	educ <- c(1,5,7)
	
	ret <- c(1,3,5)
	
		match2 <- matchit(treated ~factor(retireW1)+factor(retireW2)+factor(retireW3)+
	                   pressure+sex+unemployed+factor(income)+factor(education)+age,data = dat_match,
                    method = "cem", cutpoints = list(age = ag, 
                    pressure = press), 
                    grouping = list(c(retireW1 = ret.grp, retireW2 = ret.grp, retireW3 = ret.grp, 
                    education = educ.grp, income = inc.grp, sex = s, unemployed = emp)), k2k = TRUE)
	
	
	
	
	
	#summary(match)
	#plot(match)


	dta_m<-match.data(match)
	matchID<-dta_m$id[dta_m$treated==0] #store ID of matched control group
	
	dta_m2<-match.data(match2)
	matchID2<-dta_m2$id[dta_m2$treated==0] #store ID of matched control group

## calculate quantities of interest for plotting

#dat.balance$partygroup<-recode(dat.balance$partyAll,'"V"=1;else=0')
dat.balance <- dat.balance %>%
	  mutate(partygroup = case_when(partyAll == "V" ~ 1,
	                             TRUE ~ 0))
dat.balance$partygroup<-factor(dat.balance$partygroup)

#dat.balance$partygroup_matched<-recode(dat.balance$partyAll,'"V"=1;else=NA')
dat.balance <- dat.balance %>%
	  mutate(parytgroup_matched = case_when(partyAll == "V" ~ 1,
	                             TRUE ~ NA_real_))
dat.balance$partygroup_matched[dat.balance$id%in%matchID] <- 0

#dat.balance$partygroup_matched2<-recode(dat.balance$partyAll,'"V"=1;else=NA')
dat.balance <- dat.balance %>%
	  mutate(partygroup_matched2 = case_when(partyAll == "V" ~ 1,
	                            TRUE ~ NA_real_))
dat.balance$partygroup_matched2[dat.balance$id%in%matchID2] <- 0

#calculate quantities of interest
SE=function(x)sd(x)/sqrt(length(x)-1) #function to calculate naive SEs
out<-aggregate(retire_cutClean~partygroup*time,data=dat.balance,mean)
out.se<-aggregate(retire_cutClean~partygroup*time,data=dat.balance,SE)

outm<-aggregate(retire_cutClean~partygroup_matched*time,data=dat.balance,mean)
outm.se<-aggregate(retire_cutClean~partygroup_matched*time,data=dat.balance,SE)

outm2<-aggregate(retire_cutClean~partygroup_matched2*time,data=dat.balance,mean)
outm2.se<-aggregate(retire_cutClean~partygroup_matched2*time,data=dat.balance,SE)


##########
# Plotting matches tables

library(tidyverse)

retire_psm <- match.data(object = match) %>%
  filter(subclass %in% c("22", "20", "115")) %>%
  select(-time, -treated, -weights) %>%
  mutate(pressure = round(pressure, digits = 2),
         distance = round(distance, digits = 2)) %>%
  rename(group = subclass,
         party = partyAll,
         educ = education,
         unemp = unemployed,
         dist = distance,
         pres = pressure,
         inc = income,
         retW1 = retireW1,
         retW2 = retireW2, 
         retW3 = retireW3) %>%
  arrange(group) %>% 
  select(group, id, dist, pres, educ, sex, unemp, age, inc, retW1, retW2, retW3, party)

table_3 <- knitr::kable(retire_psm, caption = "PSM Match on Retirement Policy Responses")

retire_cem <- match.data(object = match2) %>%
  filter(subclass %in% c("2", "17", "9")) %>%
  select(-time, -treated, -weights) %>%
  mutate(pressure = round(pressure, digits = 2)) %>%
  rename(group = subclass,
         party = partyAll,
         educ = education,
         unemp = unemployed,
         pres = pressure,
         inc = income,
         retW1 = retireW1,
         retW2 = retireW2, 
         retW3 = retireW3) %>%
  arrange(group) %>% 
  select(group, id, pres, educ, sex, unemp, age, inc, retW1, retW2, retW3, party)

table_4 <- knitr::kable(retire_cem, caption = "CEM Match on Retirement Policy Responses")
  
```

```{r}
table_3
```

```{r}
table_4
```

```{r, include=FALSE}

## PROPENSITY IMBALANCE


plot(summary(match))
```

```{r, include=FALSE}

## CEM IMBALANCE

plot(summary(match2))
```

```{r, include=FALSE}

## UNMATCHED IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

#dat_bal <- dat.balance
#dat_bal$treated<-recode(dat_bal$partyAll,'"V"=1;"DF"=1;else=0')


cov_means <- dat_match %>%
  group_by(treated) %>%
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var <- dat_match %>%
  group_by(treated) %>%
   dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre0 <- cov_means %>% 
  left_join(cov_var, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Before matching')


ggplot(bal_pre0, aes(x = var, y = SMD)) + 
  geom_bar(stat = 'identity',
           col = 'black') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()


```

```{r, include=FALSE}

## PROPENSITY SCORE IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

cov_means <- dta_m %>%
  group_by(treated) %>%
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var <- dta_m %>%
  group_by(treated) %>%
   dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre <- cov_means %>% 
  left_join(cov_var, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Propensity Score')


plot <- ggplot(bal_pre, aes(x = var, y = SMD)) + 
  geom_bar(stat = 'identity',
           col = 'black') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()

# plot
```

```{r, include=FALSE}

## PROPENSITY SCORE VS CEM IMBALANCE

library(tidyverse)
library(broom)
library(StatMatch)

#dat_bal <- dat.balance
#dat_bal$treated<-recode(dat_bal$partyAll,'"V"=1;"DF"=1;else=0')


cov_means2 <- dta_m2 %>%
  group_by(treated) %>%
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var2 <- dta_m2 %>%
  group_by(treated) %>%
   dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)


## Combine 
## Calculate standardized mean difference (SMD)

bal_pre2 <- cov_means2 %>% 
  left_join(cov_var2, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'CEM')


#total_bal <- left_join(bal_pre,bal_pre2,by = c('var', 'type'))

total_bal <- rbind(bal_pre,bal_pre2)



plot <- ggplot(total_bal, aes(x = var, y = SMD, fill = type)) + 
  geom_bar(stat = 'identity',
           col = 'black', position = 'dodge') +
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw()

# plot

```

To ensure that the CEM matches were all closely matched beyond the ones in table 4, figure 3 shows the normalized difference of all the variables used before matching, using CEM matching, and using PSM matching. The CEM matching clearly outperformed the PSM and no matching, reducing the difference in most variables to 0. Pressure and age were the only variables that had shifts in differences.

```{r, warning=FALSE, fig.cap = "Comparison of variable imbalance among the respondents' to retirement policy. CEM matching results in a better match than the PSM matching used in the original paper."}

# Figure 3

## COMPLETE COMPARISON

library(StatMatch)
## Step 1: MAHALANOBIS DISTANCE MATCHING 
## Covariate Data Matrix 

big_X <- dat_match %>% 
    dplyr::select(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Only treated units 

X_treated <- dat_match %>%
  filter(treated == 1) %>%
    dplyr::select(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Only control units 

X_control <- dat_match %>%
  filter(treated == 0) %>%
  dplyr::select(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education"))) %>% 
  as.matrix()

## Calculate Mahalanobis Distance between all observation-pairs 
## using mahalanobis.dist function from the StatMatch package 

M <- mahalanobis.dist(data.x = X_treated, 
                      data.y = X_control,
                      vc = var(big_X))

## For each treated unit, find closest control unit 
## in terms of Mahalanobis Distance 

match_ids <- unlist(apply(M, 1, function(i){which(i == min(i))[1]}))

## Create matched dataset 

matched_data <- rbind(subset(dat_match, treated == 1),
                      subset(dat_match, treated == 0)[match_ids,])

matched_data_cem <- dta_m2

matched_data_prop <- dta_m


## Check balance after matching ! 

cov_means_matched <- matched_data %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched <- matched_data %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post <- cov_means_matched %>% 
  left_join(cov_var_matched, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'Mahalanobis')
#####
cov_means_matched2 <- matched_data_cem %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched2 <- matched_data_cem %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post2 <- cov_means_matched2 %>% 
  left_join(cov_var_matched2, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'CEM')

####
cov_means_matched3 <- matched_data_prop %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(mean_control = 2,
         mean_treated = 3)


cov_var_matched3 <- matched_data_prop %>% 
  group_by(treated) %>% 
  dplyr::summarise(across(c("education", "pressure", "unemployed", "sex", "age", "income",factor("retireW1"), factor("retireW2"), factor("retireW3"), factor("income"), factor("education")), mean), .groups = "drop") %>%
  t() %>%
  as_tibble(rownames = 'var') %>%
  slice(-1) %>%
  rename(var_control = 2,
         var_treated = 3)

bal_post3 <- cov_means_matched3 %>% 
  left_join(cov_var_matched3, by = "var") %>%
  mutate(mean_diff = mean_treated - mean_control,
         SMD = mean_diff / sqrt((var_control + var_treated)/2)) %>%
  mutate(type = 'PSM')

## Compare balance before and after matching visually 

ggplot(rbind(bal_pre0, bal_post2, bal_post3), aes(x = var, y = SMD, 
                                     col = type, group = type, shape = type)) +
  geom_point() + 
  coord_flip() + 
  labs(y = 'Normalized difference',
       x = 'Variable') + 
  theme_bw() + 
  geom_hline(yintercept = 0,
             linetype = 'dotted') + 
  theme(legend.position = 'bottom',
        legend.title = element_blank())
```

\newpage

Figure 4 shows the original policy support response averages by the PSM matched, the control, and the treated groups. We add our CEM matched averages to show potentially contradicting results to original findings. The more strictly matched pairs using our cut-points showed that the average policy support of the matched respondents might actually be much higher at about 0.80. The treated Liberal respondents have a slightly lower support of abolishing early retirement between 0.65 and 0.70. The average of the treated respondents only increase their support of the policy to the 0.8 after their party shifts positions on the policy. This suggests that the original support for abolishing early retirement among similar respondents to the treatment was much higher than originally thought to be. 

```{r, fig.cap = "Respondents' support for the policy to abolish early retirement. Vertical dotted line represents the date in which the treatment party shifted their policy support. Confidence Inervals for each point are at the 95% level."}
# Figure 4

##########
#PLOTTING

c1="grey70"
c2 = "dodgerblue"
c3 = "grey50"
#pdf('retire_plot.pdf',width=6,height=4)

  par(omd=c(0.05,1,0,1),mai=c(1,1,.1,.1),mfrow=c(1,1))
  y=out[out$partygroup==1,3]
  se=out.se[out.se$partygroup==1,3]
  plot(x,y,ylim=c(.4,1),xlim=c(x[1],as.Date("2011-08-01")),main="",xaxt="n",yaxt="n",ylab="",xlab="")
  segments(x,y+1.96*se,x,y-1.96*se)
  lines(x,y,lty=1)
  points(x,y,pch=21,bg="black")
  y=out[out$partygroup==0,3]
  se=out.se[out.se$partygroup==0,3]
  lines(x,y,lty=1,col=c1)
  points(x,y,pch=24,bg=c1,col=c1)
  segments(x,y+1.96*se,x,y-1.96*se,col=c1)
  y=outm[outm$partygroup==0,3]
  se=outm.se[outm.se$partygroup==0,3]
  lines(x,y,lty=1,col=c3)
  points(x,y,pch=21,bg=c3,col=c3)
  segments(x,y+1.96*se,x,y-1.96*se,col=c3)
  y=outm2[outm2$partygroup==0,3]
  se=outm2.se[outm2.se$partygroup==0,3]
  lines(x,y,lty=1,col=c2)
  points(x,y,pch=21,bg=c2,col=c2)
  segments(x,y+1.96*se,x,y-1.96*se,col=c2)

#redraw CI's of treatment group (for overlap)
  y=out[out$partygroup==1,3]
  se=out.se[out.se$partygroup==1,3]
  segments(x,y+1.96*se,x,y-1.96*se)
  points(x,y,pch=21,bg="black")

  axis(2,las=2)
  axis(2,at=seq(0,1,.025),tck=-.015,labels=NA)
  mtext(2,text="Policy Support",line=3)
  axis.Date(1,at=seq(as.Date("2010-03-01"),as.Date("2011-08-01"), by="month"),format="%b",las=1,cex.axis=1)
  axis.Date(1,at=c(as.Date("2010-03-01"),as.Date("2011-01-01")),
            format="%Y",las=1,cex.axis=1,line=1,tick="FALSE")
  abline(v=as.Date("2011-01-01"),lty=2)
  legend('bottomright',pch=c(21,24,21),pt.bg=c("black",c1,c3,c2),
         col=c("black",c1,c3,c2),legend=c("Treated","Control (all)","Matched(PSM)","Matched (CEM)"),cex=.8,pt.cex=1,bty='n')

#dev.off()

```

# Discussion

Next, we reviewed evidence for reverse causation. Increased public support for these policies could have plausibly led to parties changing platforms in the first place. Although the initial level of support for retirement policies was higher than expected among the CEM-matched control group, which indicates that public opinion was higher for non-liberal partisans prior to the liberal parties’ policy changes, it is inconclusive whether this increased public support influenced liberal parties to target retirement policy specifically. Although panel data provides a unique opportunity to control for reverse causation through lagged independent variables, the short timeframe of the panel data disables this analysis.

# Conclusion

In this paper, we have replicated and analyzed Slothuus and Bisgaard’s study of the impact of party policy changes towards unemployment benefits and early retirement on public support. Through further examination of their use of propensity score matching, we found that the resulting matches had notable differences between values. This imbalance can produce biased results and therefore should be improved. Using coarsened exact matching and strategic cut points based on distribution analysis, we were able to produce matches with normalized differences between pairs of matched observations that were consistently smaller than those from propensity score matching. Using this improvement, we found that the original trends of public support were not as consistent as the original paper demonstrates. When replicating the figures from the original paper, we found that the policy for unemployment benefits corroborated the analysis of the paper, while the policy for early retirement did not. With further analysis of the data, we found that there exists a large imbalance between unemployed and employed observations. Therefore, the trend of policy support regarding unemployment benefits represents the beliefs of only the employed public. The early retirement public support trends found using coarsened matching showed a higher approval of the policy even before the political parties shift in policies. With the improvement of matching and the discovery of a biased dataset, we find that political parties’ influence on shifting public support has a lesser impact than previously shown. 

\newpage

# Bibliography

